Spatial Multivariate Trees for Big Data Bayesian Regression

High resolution geospatial data are challenging because standard geostatistical models based on Gaussian processes are known to not scale to large data sizes. While progress has been made towards methods that can be computed more efficiently, considerably less attention has been devoted to methods for large scale data that allow the description of complex relationships between several outcomes recorded at high resolutions by different sensors. Our Bayesian multivariate regression models based on spatial multivariate trees (SpamTrees) achieve scalability via conditional independence assumptions on latent random effects following a treed directed acyclic graph. Information-theoretic arguments and considerations on computational efficiency guide the construction of the tree and the related efficient sampling algorithms in imbalanced multivariate settings. In addition to simulated data examples, we illustrate SpamTrees using a large climate data set which combines satellite data with land-based station data. Software and source code are available on CRAN at https://CRAN.R-project.org/package=spamtree.


Introduction
It is increasingly common in the natural and social sciences to amass large quantities of geo-referenced data. Researchers seek to use these data to understand phenomena and make predictions via interpretable models that quantify uncertainty taking into account the spatial and temporal dimensions. Gaussian processes (GP) are flexible tools that can be used to characterize spatial and temporal variability and quantify uncertainty, and considerable attention has been devoted to developing GP-based methods that overcome their notoriously poor scalability to large data. The literature on scaling GPs to large scale is now extensive. We mention low-rank methods (Quiñonero-Candela and Rasmussen, 2005;Snelson and Ghahramani, 2007;Banerjee et al., 2008;Cressie and Johannesson, 2008); their extensions (Low et al., 2015;Ambikasaran et al., 2016;Huang and Sun, 2018;Geoga et al., 2020); methods that exploit special structure or simplify the representation of multidimensional inputs-for instance, a Toeplitz structure of the covariance matrix scales GPs to big time series data, and tensor products of scalable univariate kernels can be used for multidimensional inputs (Gilboa et al., 2015;Moran and Wheeler, 2020;Loper et al., 2020;Wu et al., 2021). These methods may be unavailable or perform poorly in geostatistical settings, which focus on small-dimensional inputs, i.e. the spatial coordinates plus time. In these scenarios, low-rank methods oversmooth the spatial surface (Banerjee et al., 2010), Toeplitz-like structures are typically absent, and so-called separable covariance functions obtained via tensor products poorly characterize spatial and temporal dependence. To overcome these hurdles, one can use covariance tapering and domain partitioning (Furrer et al., 2006;Kaufman et al., 2008;Sang and Huang, 2012;Stein, 2014;Katzfuss, 2017) or composite likelihood methods and sparse precison matrix approximations (Vecchia, 1988;Rue and Held, 2005;Eidsvik et al., 2014); refer to Sun et al. (2011), Banerjee (2017), Heaton et al. (2019) for reviews of scalable geostatistical methods.
Additional difficulties arise in multivariate (or multi-output) regression settings. Multivariate geostatistical data are commonly misaligned, i.e. observed at non-overlapping spatial locations (Gelfand et al., 2010). Figure 1 shows several variables measured at nonoverlapping locations, with one measurement grid considerably sparser than the others. In these settings, replacing a multi-output regression with separate single-output models is a valid option for predicting outcomes at new locations. While single-output models may some-times perform equally well or even outperform multi-output models, they fail to characterize and estimate cross-dependences across outputs; testing the existence of such dependences may be scientifically more impactful than making predictions. This issue can be solved by modeling the outputs via latent spatial random effects thought of as a realization of an underlying multivariate GP and embedded in a larger hierarchical model.
Unfortunately, GP approximations that do not correspond to a valid stochastic process may inaccurately characterize uncertainty, as the models used for estimation and interpolation may not coincide. Rather than seeking approximations to the full GP, one can develop valid standalone spatial processes by introducing conditional independence across spatial locations as prescribed by a sparse directed acyclic graph (DAG). These models are advantageous because they lead to scalability by construction; in other words, posterior computing algorithms for these methods can be interpreted not only as approximate algorithms for the full GP, but also as exact algorithms for the standalone process. domain partitioning; in this case, too, efficiencies follow from the properties of the chosen DAG. There is a rich literature on Gaussian processes and recursive partitioning, see e.g Ferreira and Lee (2007); Gramacy and Lee (2008); Fox and Dunson (2012); in geospatial contexts, in addition to the GMRF-based method of Nychka et al. (2015), multi-resolution approximations (MRA; Katzfuss, 2017) replace an orthogonal basis decomposition with approximations based on tapering or domain partitioning and also have a DAG interpretation (Katzfuss and Guinness, 2021).
Considerably less attention has been devoted to process-based methods that ensure scalability in multivariate contexts, with the goal of modeling the spatial and/or temporal variability of several variables jointly via flexible cross-covariance functions (Genton and Kleiber, 2015). When scalability of GP methods is achieved via reductions in the conditioning sets, including more distant locations is thought to aid in the estimation of unknown covariance parameters (Stein et al., 2004). However, the size of such sets may need to be reduced excessively when outcomes are not of very small dimension. One could restrict spatial coverage of the conditioning sets, but this works best when data are not misaligned, in which case all conditioning sets will include outcomes from all margins; this cannot be achieved for misaligned data, leading to pathological behavior. Alternatively, one can model the multivariate outcomes themselves as a DAG; however this may only work on a case-by-case basis. Similarly, recursive domain partitioning strategies work best for data that are measured uniformly in space as this guarantees similarly sized conditioning sets; on the contrary, recursive partitioning struggles in predicting the outcomes at large unobserved areas as they tend to be associated to the small conditioning sets making up the coarser scales or resolutions.
In this article, we solve these issues by introducing a Bayesian regression model that encodes spatial dependence as a latent spatial multivariate tree (SPAMTREE); conditional independence relations at the reference locations are governed by the branches in a treed DAG, whereas a map is used to assign all non-reference locations to leaf nodes of the same DAG. This assignment map controls the nature and the size of the conditioning sets at all locations; when severe restrictions on the reference set of locations become necessary due to data size, this map is used to improve estimation and predictions and overcome common issues in standard nearest-neighbor and recursive partition methods while maintaining the desirable recursive properties of treed DAGs. Unlike methods based on defining conditioning sets based solely on spatial proximity, SPAMTREES scale to large data sets without excessive reduction of the conditioning sets. Furthermore, SPAMTREES are less restrictive than methods based on recursive partitioning and can be built to guarantee similarly-sized conditioning sets at all locations.
The present work adds to the growing literature on spatial processes defined on DAGs by developing a method that targets efficient computations of Bayesian multivariate spatial regression models. SPAMTREES share similarities with MRAs (Katzfuss, 2017); however, while MRAs are defined as a basis function expansion, they can be represented by a treed graph of a SPAMTREE with full "depth" as defined later (the DAG on the right of Figure 2), in univariate settings, and "response" models. All these restrictions are relaxed in this article. In considering spatial proximity to add "leaves" to our treed graph, our methodology also borrows from nearest-neighbor methods (Datta et al., 2016a). However, while we use spatial neighbors to populate the conditioning sets for non-reference locations, the same cannot be said about reference locations for which the treed graph is used instead. Our construction of the SPAMTREE process also borrows from MGPs on tessellated domains (Peruzzi et al., 2020); however, the treed DAG we consider here induces markedly different properties on the resulting spatial process owing to its recursive nature. Finally, a contribution of this article is in developing self-contained sampling algorithms which, based on the graphical model representation of the model, will not require any external libraries.
The article builds SPAMTREES as a standalone process based on a DAG representation in Section 2. A Gaussian base process is considered in Section 3 and the resulting properties outlined, along with sampling algorithms. Simulated data and real-world applications are in Section 4; we conclude with a discussion in Section 5. The Appendix provides more in-depth treatment of several topics and additional algorithms.

Spatial Multivariate Trees
Consider a spatial or spatiotemporal domain D. With the temporal dimension, we have D ⊂ ℜ d × 0, ∞ , otherwise D ⊂ ℜ d . A q-variate spatial process is defined as an uncountable set of random variables w ℓ : ℓ ∈ D , where w ℓ is a q × 1 random vector with elements w i ℓ for i = 1, 2,…, q, paired with a probability law P defining the joint distribution of any finite sample from that set. Let ℓ 1 , ℓ 2 , …, ℓ n ℒ = ℒ ⊂ D be of size n ℒ . The n ℒ q × 1 random vector w ℒ = w ℓ 1 ⊤ , w ℓ 2 ⊤ , …w ℓ n ℒ ⊤ ⊤ has joint density p w ℒ . After choosing an arbitrary order of the locations, p w ℒ = i = 1 n ℒ p w ℓ i | w ℓ 1 , …, w ℓ i − 1 , where the conditioning set for each w ℓ i can be interpreted as the set of nodes that have a directed edge towards w ℓ i in a DAG. Some scalable spatial processes result from reductions in size of the conditioning sets, following one of several proposed strategies (Vecchia, 1988;Stein et al., 2004;Gramacy and Apley, 2015;Datta et al., 2016a;Katzfuss and Guinness, 2021;Peruzzi et al., 2020). Accordingly, where Pa[ℓ i ] is the set of spatial locations that correspond to directed edges pointing to ℓ i in the DAG. If Pa[ℓ i ] is of size J or less for all i = 1, …, n ℒ , then w Pa ℓ i is of size J q . Methods that rely on reducing the size of parent sets are thus negatively impacted by the dimension q of the multivariate outcome; if q is not very small, reducing the number of parent locations J may be insufficient for scalable computations. As an example, an NNGP model has Pa[ℓ i ] = N (ℓ i ), where N (·) maps a location in the spatial domain to its neighbor set. It is customary in practice to consider J q = m ≤ 20 for accurate and scalable estimation and predictions We represent the ith component of the q × 1 vector w ℓ as w ℓ, ξ i , where ξ i = ξ i1 , …, ξ ik ⊤ ∈ Ξ for some k and Ξ serves as the k-dimensional latent spatial domain of variables. The q-variate process w ℓ is thus recast as w ℓ, ξ : ℓ, ξ ∈ D × Ξ , with ξ representing the latent location in the domain of variables. We can then write (1) as where ℒ * = ℓ i * i = 1 n ℓ * , ℓ i * ∈ D × Ξ = D * , and w(·) is a univariate process on the expanded domain D * . This representation is useful as it provides a clearer accounting of the assumed conditional independence structure of the process in a multivariate context.

Constructing Spatial Multivariate DAGs
We now introduce the necessary terminology and notation, which are the basis for later detailing of estimation and prediction algorithms involving SPAMTREES. The specifics for building treed DAGs with user-specified depth are in Section 2.1.1, whereas Section 2.1.2 gives details on cherry picking and its use when outcomes are imbalanced and misaligned.
The three key components to build a SPAMTREE are (i) a treed DAG G with branches and leaves on M levels and with depth δ ≤ M; (ii) a reference set of locations S; (iii) a cherry picking map. The graph is G = V , E where the nodes are V = v 1 , …, v m V = A ∪ B, A ∩ B = 0. We separate the nodes into reference A and non-reference B nodes, as this will aid in showing that SPAMTREES lead to standalone spatial processes in Section 2.2. The reference or branch nodes are A = a 1 , …, a m A = A 0 ∪ A 1 ∪ ⋯ ∪ A M − 1 , where A i = a i, 1 , …, a i, m i for all i = 0,…, M − 1 and with A i ∩ A j = 0 if i ≠ j. The non-reference or leaf nodes are B = b 1 , …, b m B , A ∩ B = 0. We also denote V r = A r for r = 0,…, M -1 and V M = B. The edges are E = Pa v ⊂ V : v ∈ V and similarly Ch v = v′ ∈ V : v ∈ Pa v′ . The reference set S is partitioned in M levels starting from zero, and each level is itself partitioned into reference subsets:S = ∪ r = 0 M − 1 S r = ∪ r = 0 M − 1 ∪ i = 1 m i S ri , where S ri ∩ S r′i′ = 0 if r ≠ r′ or i ≠ i′ and its complement set of non-reference or other locations U = D * \S. The cherry picking map is η: D * V and assigns a node (and therefore all the edges directed to it in G) to any location in the domain, following a user-specified criterion.
2.1.1. BRANCHES AND LEAVES-For a given M and a depth δ ≤ M , we impose a treed structure on G by assuming that if v ∈ A i and i > M − δ = M δ then there exists a sequence of nodes v r M δ , …, v r i − 1 such that v r j ∈ A j for j = M δ , …, i − 1 and is the tree root and is such that Pa v 0 = 0 for all v 0 ∈ A 0 . The depth δ determines the number of levels of G (from the top) across which the parent sets are nested. Choosing δ = 1 implies that all nodes have a single parent; choosing δ = M implies fully nested parent sets (i.e. if v i ∈ Pa v j then Pa v i ⊂ Pa v j for all v i , v i ∈ V . The m i elements of A i are the branches at level i of G and they have i -M δ parents if the current level i is above the depth level M δ and 1 parent otherwise. We refer to terminal branches as nodes v ∈ A such that Ch v ⊂ B. For all choices of δ, v ∈ A i , v′ ∈ A j and v ∈ Pa v′ implies i < j; this guarantees acyclicity.
As for the leaves, for all v ∈ B we assume Pa v = v r M δ , …, v r k for some integer sequence r M δ , …, r k and v r i ∈ A i with i ≥ M δ . We allow the existence of multiple leaves with the same parent set, i.e. there can be k and b i 1 , …, b i k such that for all i 2 ,…, i k ,Pa b i ℎ = Pa b i 1 .
Acyclicity of G is maintained as leaves are assumed to have no children. Figure 2 represents the graph associated to SPAMTREES with different depths.

CHERRY PICKING VIA η(·)-
The link between G, S and U is established via the map η: D * V which associates a node in G to any location ℓ * in the expanded domain D * : This is a many-to-one map; note however that all locations in S ij are mapped to a ij : by calling η X = η ℓ * : ℓ * ∈ X then for any i = 0,…, M -1 and any j = 1,…, m i we have η S ij = η A S ij = a ij . SPAMTREES introduce flexibility by cherry picking the leaves, i.e. using η B : U B, the restriction of η to U. Since each leaf node b j determines a unique path in G ending in b j , we use η B to assign a convenient parent set to w(u), u ∈ U, following some criterion.
For example, suppose that u = ℓ, ξ s meaning that w u = w ℓ, ξ s is the realization of the sth variable at the spatial location ℓ, and we wish to ensure that Pa[w(u)] includes realizations of the same variable. Denote T = v ∈ A: Ch v ⊂ B as the set of terminal branches of G. Then we find ℓ, ξ s opt = argmin ℓ′, ξ′ = ξ s ∈ η A −1 T d ℓ′, ℓ where d(·,·) is the Euclidean distance. Since ℓ, ξ s opt ∈ S ij for some i, j we have η A ℓ, ξ s opt = a ij . We then set η B (u) = b k where Pa[b k ] = {a ij }. In a sense a ij is the terminal node nearest to u; having defined η B in such a way forces the parent set of any location to include at least one realization of the process from the same variable. There is no penalty in using D * = D × Ξ as we can write p w u | Pa w u = p w ℓ, ξ 1 , …, ℓ, ξ q | Pa w u = s = 1 q p w ℓ, ξ s | w ℓ, ξ 1 , …, w ℓ, ξ s − 1 , Pa w ℓ , which also implies that the size of the parent set may depend on the variable index. Assumptions of conditional independence across variables can be encoded similarly. Also note that any specific choice of η B induces a partition on U; let U j = u ∈ U: η B u = b j , then clearly U = ∪ j = 1 m U U j with U i ∩ U j = 0 inactive, resulting in m U = m B . We will thus henceforth assume that m U = m B without loss of generality.

SPAMTREES as a Standalone Spatial Process
We define a valid joint density for any finite set of locations in D * satisfying the Kolmogorov consistency conditions in order to define a valid process. We approach this problem analogously to Datta et al. (2016a) and Peruzzi et al. (2020). Enumerate each of the m S reference subsets as S i = s i 1 , …, s i n i where i 1 , …, i n i ⊂ 1, …, n S , and each of the m U non-reference subsets as U i = u i 1 , …, u i n i where i 1 , …, i n i ⊂ 1, …, n U . Then Then take w i = w ℓ i 1 , …, w ℓ i n i ⊤ as the n i × 1 random vector with elements of w ℓ for each ℓ ∈ V i . Denote w i = w η −1 Pa v i . Then which is a proper multivariate joint density since G is acyclic (Lauritzen, 1996). All locations inside U j always share the same parent set, but a parent set is not necessarily unique to a single U j . This includes as a special case a scenario in which one can assume in this case each location corresponds to a leaf node. To conclude the construction, for any finite subset of spatial locations ℒ ⊂ D we can let U = ℒ\S and obtain p w ℒ = p w U w S p w S s i ∈ S\ℓ d w s i , leading to a well-defined process satisfying the Kolmogorov conditions (see Appendix A).

POSITIONING OF SPATIAL LOCATIONS IN CONDITIONING SETS-
In spatial models based on sparse DAGs, larger conditioning sets yield processes that are closer to the base process p in terms of Kullback-Leibler divergence (Banerjee, 2020;Peruzzi et al., 2020), denoted as KL p | | ⋅ . The same results cannot be applied directly to SPAMTREES given the treed structure of the DAG. For a given S, we consider the distinct but related issues of placing individual locations into reference subsets (1) at different levels of the treed hierarchy; (2) within the same level of the hierarchy.
The proof proceeds by an "information never hurts" argument (Cover and Thomas, 1991).
Intuitively, this shows that there is a penalty associated to positioning reference locations at higher levels of the treed hierarchy. Increasing the size of the reference set at the root augments the conditioning sets at all its children; since this is not true when the increase is at a branch level, the KL divergence of p 0 from p is smaller than the divergence of p 1 from the same density. In other words there is a cost of branching in G which must be justified by arguments related to computational efficiency. The above proposition also suggests populating near-root branches with locations of sparsely-observed outcomes. Not doing so in highly imbalanced settings may result in possibly too restrictive spatial conditional independence assumptions.
While we do not target the estimation of these quantities, this result is helpful in designing SPAMTREES as it suggests placing a new reference location s * in the reference subset least uncertain about the realization of the process at s * . We interpret this as justifying recursive domain partitioning on S in spatial contexts in which local spatial clusters of locations are likely less uncertain about process realization in the same spatial region. In the remainder of this article, we will consider a given reference set S which typically will be based on a subset of observed locations; the combinatorial problem of selecting an optimal S (in some sense) is beyond the scope of this article. If S is not partitioned, it can be considered as a set of knots or "sensors" and one can refer to a large literature on experimental design and optimal sensor placement (see e.g. Krause et al., 2008, and references therein). It might be possible to extend previous work on adaptive knot placement (Guhaniyogi et al., 2011), but this will come at a steep cost in terms of computational performance.

Bayesian Spatial Regressions Using SPAMTREES
Suppose we observe an l-variate outcome at spatial locations ℓ ∈ D ⊂ ℜ d which we wish to model using a spatially-varying regression model: where y j ℓ is the j-th point-referenced outcome at ℓ, x j ℓ is a p j × 1 vector of spatially referenced predictors linked to constant coefficients β j , ε j ℓ iid N 0, τ j 2 is the measurement error for outcome j, and z jk ℓ is the k-th (of q) covariates for the j-th outcome modeled with spatially-varying coefficient w ℓ, ξ k , ℓ ∈ D,ξ k ∈ Ξ. This coefficient w ℓ, ξ k corresponds to the k-th margin of a q-variate Gaussian process w ℓ : ℓ ∈ D denoted as w ℓ GP 0, C θ ⋅ , ⋅ with cross-covariance C θ indexed by unknown parameters θ which we omit in notation for simplicity. A valid cross-covariance function is defined as C θ : D × D ℳ q × q , where ℳ q × q is a subset of the space of all q × q real matrices ℜ q × q . It must satisfy C ℓ, ℓ′ = C ℓ′, ℓ ⊤ for any two locations ℓ, ℓ′ ∈ D and i = 1 n j = 1 n We replace the full GP with a Gaussian SPAMTREE for scalable computation considering the q-variate multivariate Gaussian process w(·) as the base process. Since the (i, j)-th entry of C ℓ, ℓ′ is C ℓ, ℓ′ i, j = Cov w i ℓ , w j ℓ′ , i.e. the covariance between the i-th and j-th elements of w ℓ at ℓ and ℓ′, we can obtain a covariance function on the augmented domain C * : D* × D* ℜ as C * ℓ, ξ , ℓ′, ξ′ = C ℓ, ℓ′ i, i′ where ξ and ξ′ are the locations in Ξ of variables i and j, respectively. Apanasovich and Genton (2010) use a similar representation to build valid cross-covariances based on existing univariate covariance functions; their approach amounts to considering ξ or ξ − ξ′ as a parameter to be estimated. Our approach can be based on any valid cross-covariance as we may just set Ξ = 1,…, q . Refer to e.g. Genton and Kleiber (2015) for an extensive review of cross-covariance functions for multivariate processes. Moving forward, we will not distinguish between C * and C.
The linear multivariate spatially-varying regression model (7) allows the l outcomes to be observed at different locations; we later consider the case l = q and Z ℓ = I q resulting in a multivariate space-varying intercept model.

Gaussian SPAMTREES
Enumerate the set of nodes as V = v 1 , …, v m V , m V = m S + m U and denote w i = w η −1 v i , C ij as the n i × n j covariance matrix between w i and w j , C i , [i] the n i × J i covariance matrix between w i and w [i] , C i the n i × n i covariance matrix between w i and itself, and C [i] the J i × J i covariance matrix between w [i] and itself. A base Gaussian process induces implying that the joint density p w S is multivariate normal with covariance C S and precision matrix C S −1 . At U we have p w U | w S = j: v j ∈ B N w j | H j w j , R j , where H j and R j are as in (8). All quantities can be computed using the base cross-covariance function. Given that the p densities are Gaussian, so will be the finite dimensional distributions.
The treed graph G leads to properties which we analyze in more detail in Appendix B and summarize here. For two nodes v i , v j ∈ V denote the common descendants as If v i ∈ Pa v j denote H i→j and H \i j as the matrix obtained by subsetting H j to columns corresponding to v i , or to Pa v j \ v j , respectively. Similarly define w i j = w i and w \i j . As a special case, if the tree depth is

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript where I ij is the (i, j) block of an identity matrix with n S + n U rows and is nonzero if and only if i = j. We thus obtain that the number of nonzero elements of C −1 is If δ > 1, the size of C [i] is larger for nodes v i at levels of the treed hierarchy farther proceeds more cheaply by recursively applying the following:

INDUCED COVARIANCE-Define
and v i ℎ ∈ Pa v i ℎ + 1 . The longest path P k j is such that if v k ∈ A r k and v j ∈ A r j then P k j = r j − r k + 1. The shortest path P k j is the path from v k to v j with minimum number of steps. We denote the longest path from the root to v j as P 0 j ; this corresponds to the full set of ancestors of v j , and Pa v j ⊂ P 0 j . For two nodes v i and v j we have Pa v i ∩ Pa v j ⊂ P 0 i ∩ P 0 j . We define the concestor between v i and v j as con v i , v j = argmax v k ∈ V k: P k i ∩ P k j ≠ 0 i.e. the last common ancestor of the two nodes.
Take the path P M δ j in G from a node at A M δ leading to v j . After defining the cross- Peruzzi where for s > i M δ the e s are independent zero-mean GPs with covariance K s ℓ, ℓ′ and we set K i M δ ℓ, ℓ′ = C ℓ, ℓ′ and e i M δ = w i M δ N 0, C i M δ . Take two locations ℓ, ℓ ꞌ such that v i = η ℓ , v j = η ℓ′ and let v z = con v i , v j ; if Pa v i ∩ Pa v j ≠ 0 then the above leads to where K z ℓ, ℓ′ = C ℓ, ℓ′ . If Pa v i ∩ Pa v j = 0 take the shortest paths P z i = i 1 , …, i r i and P z j = j 1 , …, j r j ; setting In particular if δ = M then Pa v i ∩ Pa v j ≠ 0 for all i, j and only (13) is used, whereas if δ = 1 then the only scenario in which (13) holds is v z = Pa v i ∩ Pa v j in which case the two are equivalent. In univariate settings, the special case in which δ = M , and hence M δ = 0, leads to an interpretation of (12) as a basis function decomposition; considering all leaf paths P j for v j ∈ B, this leads to an MRA (Katzfuss, 2017;Katzfuss and Gong, 2019). On the other hand, keeping other parameters constant, δ < M and in particular δ = 1 may be associated to savings in computing cost, leading to a trade-off between graph complexity and size of reference subsets; see Appendix B.5.

BLOCK-SPARSE CHOLESKY DECOMPOSITIONS-In recent work Jurek and Katzfuss
(2020) consider sparse Cholesky decompositions of co-variance and precision matrices for treed graphs corresponding to the case δ = M above in the context of space-time filtering; their methods involve sparse Cholesky routines on reverse orderings of C −1 at the level of individual locations. In doing so, the relationship between Cholesky decompositions and G, C −1 and the block structure in S remains somewhat hidden, and sparse Cholesky libraries are typically associated to bottlenecks in MCMC algorithms. However we note that a consequence of (9) is that it leads to a direct algorithm, for any δ, for the blockdecomposition of any symmetric positive-definite matrix Λ conforming to G, i.e. with the same block-sparse structure as C −1 . This allows us to write Λ = I − L ⊤ D I − L where I is the identity matrix, L is block lower triangular with the same block-sparsity pattern as ℋ above, and D is block diagonal symmetric positive-definite. In Appendix B.2.3 we outline Algorithm 4 which (i) makes direct use of the structure of G, (ii) computes the decomposition at blocks of reference and non-reference locations, and (iii) requires no external sparse matrix library, in particular no sparse Cholesky solvers. Along with Algorithm 5 for the block-computation of (I -L) −1 , it can be used to compute where Σ is a block-diagonal matrix; it is thus useful in computing the Gaussian integrated likelihood.

Estimation and Prediction
We introduce notation to aid in obtaining the full conditional distributions. Write (7) as where y ℓ = y j ℓ j = 1 l ⊤ The l × q matrix Z ℓ = z j ℓ ⊤ , j = 1, …, l with z j ℓ ⊤ = z jk ℓ , k = 1, …, q acts a design matrix for spatial location ℓ. Collecting all locations along the j-th margin, we build T j = ℓ 1 j , …, ℓ N j j and T = ∪ j T j . We then call y j = y j ℓ 1 j , …, y j ℓ N j j ⊤ and ε j similarly, X j = x j ℓ 1 The full observed data are y, X, Z. Denoting the number of observations as n = j = 1 l N j , Z is thus a n × qn block-diagonal matrix, and similarly w is a qn × 1 vector. We introduce the diagonal matrix D n such that diag D n = τ 1 By construction we may have η S i = v i and η S j = v j such that ℓ, ξ ∈ S i and ℓ′, ξ′ ∈ S j where ℓ′ = ℓ, ξ ≠ ξ′ and similarly for non-reference subsets. Suppose A ⊂ D × Ξ is a generic reference or non-reference subset. We denote A ⊂ D × Ξ as the set of all combinations of spatial locations of A and variables i.e. A = A D × A Ξ where A D ⊂ D is the set of unique spatial locations in A and A Ξ are the unique latent variable coordinates. By subtraction we find A − = A\A as the set of locations whose spatial location is in A but whose variable is not. Let y A = y A = y ℓ , ℓ ∈ A D ⊤ ,X A = X A = b . diag X ℓ ⊤ , ℓ ∈ A D ; values corresponding to unobserved locations will be dealt with by defining D n A as the diagonal matrix obtained from D n by replacing unobserved outcomes with zeros. Denote Z A = b . diag Z ℓ , ℓ ∈ A D and w A similarly. If A includes L unique spatial locations then y A is a L l × l vector and X A is a L l × pl matrix. In particular, Z A is a L l × Lql matrix; the subset of its columns with locations in A is denoted as Z A whereas at other locations we get Z A − . We can then separate the contribution of w A to y A from the With customary prior distributions β N 0, V β and τ j 2 Inv . Gamma a τ , b τ along with a Gaussian SPAMTREE prior on w, we obtain the posterior distribution as p w, β, τ j 2 j = 1 l , θ y ∝ p y w, β, τ j 2 j = 1 l p w θ p θ p β j = 1 l p τ j 2 .
Peruzzi We compute the full conditional distributions of unknowns in the model, save for θ; iterating sampling from each of these distributions corresponds to a Gibbs sampler which ultimately leads to samples from the posterior distribution above.
For j = 1, …, l, p τ j 2 | β, y, w = Inv . Gamma a τ, j * , b τ, j * where a τ, j * = a τ + N j /2 and b τ, j Take a node v i ∈ V . If v i ∈ A then η −1 (v i ) = S i and for v j ∈ Ch v i denote If Sampling of w at nodes at the same level r proceeds in parallel given the assumed conditional independence structure in G. It is thus essential to minimize the computational burden at levels with a small number of nodes to avoid bottlenecks. In particular computing F i c and m i c can become expensive at the root when the number of children is very large. In Algorithm 3 we show that one can efficiently sample at a near-root node v i by updating F i c and m i c via message-passing from the children of v i .

UPDATE OF θ-
The full conditional distribution of θ-which may include ξ j for j = 1,…, q or equivalently δ ij = ξ i − ξ j if the chosen cross-covariance function is defined on a latent domain of variables-is not available in closed form and sampling a posteriori can proceed Algorithm 1: Computing p(w|θ).
Algorithm 2: Sampling from the full conditional distribution of w i when δ = 1.
Algorithm 3: Sampling from the full conditional distribution of w j when δ = M.
via Metropolis-Hastings steps which involve accept/reject steps with acceptance probability α = min 1, p w | θ′ p θ′ q θ | θ′ p w | θ p θ q θ′ | θ . In our implementation, we adaptively tune the standard deviation of the proposal distribution via the robust adaptive Metropolis algorithm (RAM; Vihola, 2012). In these settings, unlike similar models based on DAG representations such as NNGPs and MGPs, direct computation via p w | θ = i N w i | H i w i , R i is inefficient as it requires computing C i −1 whose size grows along the hierarchy in G. We thus outline Algorithm 1 for computing p w | θ via (11). As an alternative we can perform the update using ratios of p y | β, θ, τ = p y | w, β, τ p w | θ dw = N y | Xβ, ZCZ ⊤ + D n using Algorithms 4 and 5 outlined in Appendix B.2.3 which require no sparse matrix library.

GRAPH COLORING FOR PARALLEL SAMPLING-An
advantage of the treed structure of G is that it leads to fixed graph coloring associated to parallel Gibbs sampling; no graph coloring algorithms are necessary (see e.g. Molloy and Reed, 2002;Lewis, 2016). Specifically, if δ = M (full depth) then there is a one to one correspondence between the M + 1 levels of G and graph colors, as evidenced by the parallel blocks in Algorithms 1 and 3. In the case δ = 1, G is associated to only two colors alternating the odd levels with the even ones. This is possible because the Markov blanket of each node at level r, with r even, only includes nodes at odd levels, and vice-versa.

COMPUTING AND STORAGE COST-
The update of τ j 2 and β can be performed at a minimal cost as typically p = j = 1 l p j is small; almost all the computation budget must be dedicated to computing p w | θ and sampling p w | y, β, τ 2 . Assume that reference locations are all observed S ⊂ T and that all reference subsets have the same size i.e. S i = N s for all i. We show in Appendix B.5 that the cost of computing SPAMTREES is O nN s 2 . As a result, SPAMTREES compare favorably to other models specifically in not scaling with the cube of the number of samples. δ does not impact the computational order, however, compared to δ = M , choosing δ = 1 lowers the cost by a factor of M or more. For a fixed reference set partition and corresponding nodes, choosing larger δ will result in stronger dependence between leaf nodes and nodes closer to the root-this typically corresponds to leaf nodes being assigned conditioning sets that span larger distances in space. The computational speedup corresponding to choosing δ = 1 can effectively be traded for a coarser partitioning of S, resulting in large conditioning sets that are more local to the leaves.

Synthetic Data
In this section we focus on bivariate outcomes (q = 2). We simulate data from model (15), setting β = 0, Z = I q and take the measurement locations on a regular grid of size 70 × 70 for a total of 4,900 spatial locations. We simulate the bivariate spatial field by sampling from the full GP using (18) as cross-covariance function; the nuggets for the two outcomes are set to τ 1 2 = 0.01 and τ 2 2 = 0.1. For j = 1, 2 we fix σ j2 = 1, α = 1, β = 1 and independently sample σ j1 U −3, 3 , ϕ j 1 ~ U(−0.1, 3), ϕ j U 0.1, 3 , ϕ U 0.1, 30 , δ 12 Exp 1 , generating a total of 500 bivariate data sets. This setup leads to empirical spatial correlations between the two outcomes smaller than 0.25, between 0.25 and 0.75, and larger than 0.75 in absolute value in 107, 330, and 63 of the 500 data sets, respectively. We introduce misalignment and make the outcomes imbalanced by replacing the first outcome with missing values at ≈50% of the spatial locations chosen uniformly at random, and then repeating this procedure for the second outcome keeping only ≈ 10% of the total locations. We also introduce almost-empty regions of the spatial domain, independently for each outcome, by replacing observations with missing values at ≈ 99% of spatial locations inside small circular areas whose center is chosen uniformly at random in [0, 1] 2 . As a result of these setup choices, each simulated data set reproduces some features of the real-world unbalanced misaligned data we consider in Section 4.2 at a smaller scale and in a controlled experiment. Figure 3 shows one of the resulting 500 data sets.
We consider SPAMTREES with δ = 1 and implement multiple variants of SPAMTREES with δ = M in order to assess their sensitivity to design parameters. Table 1 reports implementation setups and the corresponding results in all cases; if the design variable "All outcomes at ℓ" is set to "No" then a SPAMTREE is built on the D × Ξ domain. If it is set to "Yes" the DAG will be built using D only -in other words, the q margins of the latent process are never separated by the DAG if they are measured at the same spatial location. "Cherry pick same outcome" indicates whether the map η(·) should search for neighbors by first filtering for matching outcomes -refer to our discussion at Section 2.1.2. We mention here if the DAG is built using D only, then the nearest neighbor found via cherry picking will always include a realization of the same margin of w ⋅ . Finally, if SPAMTREE is implemented with "Root bias" then the reference set and the DAG are built with locations of the more sparsely observed outcome closer to root nodes of the tree as suggested by Proposition 1.  Table 1. We also include results from a non-spatial regression using Bayesian additive regression trees (BART; Chipman et al., 2010) which uses the domain coordinates as covariates in addition to a binary fixed effect corresponding to the outcome index. All methods were setup to target a compute time of approximately 15 seconds for each data set. We focused on comparing the different methods under computational constraints because (a) without constraints it would not be feasible to implement the methods for many large simulated spatial datasets; and (b) the different methods are mostly focused on providing a faster approximation to full GPs; if constraints were removed one would just be comparing the same full GP method. Table 1 reports average performance across all data sets. All Bayesian methods based on latent GPs exhibit very good coverage; in these simulated scenarios, SPAMTREES exhibit comparatively lower out-of-sample prediction errors. All SPAMTREES perform similarly, with the best out-of-sample predictive performance achieved by the SPAMTREES cherry picking based solely on spatial distance (i.e. disregarding whether or not the nearest-neighbor belongs to the same margin). Additional implementation details can be found in Appendix C.2.1. Finally, we show in Figure 4 that the relative gains of SPAMTREES compared to independent univariate NNGP model of the outcomes are increasing with the magnitude of the correlations between the two outcomes, which are only available due to the simulated nature of the data sets.

Climate Data: MODIS-TERRA and GHCN
Climate data are collected from multiple sources in large quantities; when originating from satellites and remote sensing, they are typically collected at high spatial and relatively low temporal resolution. Atmospheric and land-surface products are obtained via postprocessing of satellite imaging, and their quality is negatively impacted by cloud cover and other atmospheric disturbances. On the other hand, data from a relatively small number of land-based stations is of low spatial but high temporal resolution. An advantage of land-based stations is that they measure phenomena related to atmospheric conditions which cannot be easily measured from satellites (e.g. precipitation data, depth of snow cover).
We consider the joint analysis of five spatial outcomes collected from two sources. First, we consider Moderate Resolution Imaging Spectroradiometer (MODIS) data from the Terra satellite which is part of the NASA's Earth Observing System. Specifically, data product MOD11C3 v. 6 provides monthly Land Surface Temperature (LST) values in a 0.05 degree latitude/longitude grid (the Climate Modeling Grid or CMG). The monthly data sets cover the whole globe from 2000-02-01 and consist of daytime and nighttime LSTs, quality control assessments, in addition to emissivities and clear-sky observations. The second source of data is the Global Historical Climatology Network (GHCN) database which includes climate summaries from land surface stations across the globe subjected to common quality assurance reviews. Data are published by the National Centers of Environmental Information (NCEI) of the National Oceanic and Atmospheric Administration (NOAA) at several different temporal resolutions; daily products report five core elements (precipitation, snowfall, snow depth, maximum and minimum temperature) in addition to several other measurements.
We build our data set for analysis by focusing on the continental United States in October, 2018. The MODIS data correspond to 359,822 spatial locations. Of these, 250,874 are collected at the maximum reported quality; we consider all remaining 108,948 spatial locations as missing, and extract (1) daytime LST (LST_Day_CMG), (2) nighttime LST (LST_Night_CMG), (3) number of days with clear skies (Clear_sky_days), (4) number of nights with clear skies (Clear_sky_nights). From the GHCN database we use daily data to obtain monthly averages for precipitation (PRCP), which is available at 24,066 spatial locations corresponding to U.S. weather stations; we log-transform PRCP. The two data sources do not share measurement locations as there is no overlap between measurement locations in MODIS and GHCN, with the latter data being collected more sparsely-this is a scenario of complete spatial misalignment. From the resulting data set of size n =1,027,562 we remove all observations in a large 3 × 3 degree area in the central U.S. (from -100W to -97W and from 35N to 38N, i.e. the red area of Figure 5) to build a test set on which we calculate coverage and RMDE of the predictions.
We implement SPAMTREES using the cross-covariance function (18). Considering that PRCP is more sparsely measured and following Proposition 1, we build SPAMTREES favoring placement of GHCN locations at root nodes. We compare SPAMTREES with a Q-MGP model built on the same cross-covariance function, and two univariate models that make predictions independently for each outcome. Comparisons with other multivariate methods are difficult due to the lack of scalable software for this data size which also deals with misalignment and imbalances across outcomes. Compute times per MCMC iteration ranged from 2.4s/iteration of the multivariate Q-MGP model, to 1.5s/iteration of the univariate NNGP model. The length of the MCMC chains (30,000 for SPAMTREES and 20,000 for Q-MGP) was such that the total compute time was about the same for both models at less than 16 hours. Univariate models cannot estimate cross-covariances of multivariate outcomes and are thus associated to faster compute times; we set the length of their MCMC chains to 15,000 for a total compute time of less than 7 hours for both models. We provide additional details about the models we implemented at Appendix C. Table 2 reports predictive performance of all models, and Figure 6 maps the predictions at all locations from SPAMTREES and the corresponding posterior uncertainties. Multivariate models appear advantageous in predicting some, but not all outcomes in this real world illustration; nevertheless, SPAMTREES outperformed a Q-MGP model using the same crosscovariance function. Univariate models perform well and remain valid for predictions, but cannot estimate multivariate relationships. We report posterior summaries of θ in Appendix C.2.2. Opposite signs of σ i1 and σ j1 for pairs of variables i, j ∈ 1, …, q imply a negative relationship; however, the degree of spatial decay of these correlations is different for each pair as prescribed by the latent distances in the domain of variables δ ij . Figure 7 depicts the resulting cross-covariance function for three pairs of variables.

Discussion
In this article, we introduced SPAMTREES for Bayesian spatial multivariate regression modeling and provided algorithms for scalable estimation and prediction. SPAMTREES add significantly to the class of methods for regression in spatially-dependent data settings. We have demonstrated that SPAMTREES maintain accurate characterization of spatial dependence and scalability even in challenging settings involving multivariate data that are spatially misaligned. Such complexities create problems for competing approaches, including recent DAG-based approaches ranging from NNGPs to MGPs.
One potential concern is the need for users to choose a tree, and in particular specify the number of locations associated to each node and the multivariate composition of locations in each node. Although one can potentially estimate the tree structure based on the data, this would eliminate much of the computational speedup. We have provided theoretical guidance based on KL divergence from the full GP and computational cost associated to different tree structures. This and our computational experiments lead to practical guidelines that can be used routinely in tree building. Choosing a tree provides a useful degree of user-input to refine and improve upon an approach.
We have focused on sampling algorithms for the latent effects because they provide a general blueprint which may be used for posterior computations in non-Gaussian outcome models; efficient algorithms for non-Gaussian big geostatistical data sets are currently lacking and are the focus of ongoing research. SPAMTREES can be built on larger dimensional inputs for general applications in regression and/or classifications; such a case requires special considerations regarding domain partitioning and the construction of the tree. In particular, when time is available as a third dimension it may be challenging to build a sparse DAG with reasonable assumptions on temporal dependence. For these reasons, future research may be devoted to building sparse DAG methods combining the advantages of treed structures with e.g. Markov-type assumptions of conditional independence, and applying SPAMTREES to data with larger dimensional inputs.
p w ℒπ = p w U π w S p w S s i ∈ S\ℒ π dw s i = p w U w S p w S s i ∈ S\ℒ dw s i = p w ℒ Implying p w ℓ 1 , …, w ℓ n = p w ℓ π 1 , …w ℓ π n .
Next, take a new location location ℓ 0 ∈ D*. Call ℒ 1 = ℒ ∪ l 0 . We want to show that p w ℒ 1 dw ℓ 0 = p w ℒ .If ℓ 0 ∈ S then ℒ 1 \S = ℒ\S = U and hence p w ℓ 1 dw ℓ 0 = p w ℓ 1 \S w S p w S s i ∈ S\ℓ 1 dw s i dw ℓ 0 = p w U w S p w S s i ∈ S\ℒ dw s i = p w ℒ .
If ℓ 0 ∉ S we have p w ℒ 1 dw ℒ 0 = p w ℒ 1 \S w S p w S s i ∈ S\ℓ 1 dw s i dw ℒ 0 = p w ℒ\S ∪ ℒ 0 w S p w S s i ∈ S\ℒ dw s i dw ℒ 0 = p w ℒ 0 w ℓ\S , w S p w ℒ\S w S p w S s i ∈ S\ℒ dw s i dw ℒ 0 = p w ℒ\S w S p w S s i ∈ S\ℒ dw s i p w l 0 w S dw ℒ 0 = p w ℒ\S w S p w S s i ∈ S\ℒ dw s i = p w ℒ .

B.1 Building the precision matrix
We can represent each conditional density N w i | H i w i , R i as a linear regression on w i : where each h ij is an n i × n j coefficient matrix representing the regression of w i given w i , ω i ind N 0, R i for i = 0, 1,…, M , and each R i is an n i × n i residual covariance matrix. We where O is the matrix of zeros, whenever j ∉ j 1 , ..., j r . Using this representation, we have H i = ℎ i, j 1 , ℎ i , j 2 , ..., ℎ i , j r , which is an n i × J i block matrix formed by stacking h i,jk side by side for k = 1,…, r. Since E w i | w i = H i w i = C i, i C i −1 w i , we obtain and R i 's can be computed from the base covariance function.
In order to continue building the precision matrix, define the block matrix ℋ = ℎ ij .
We can write where ⋅ , ℎ refers to the h-th block column. More compactly using the indicator function 1 ⋅ we have h ij = 1 ∃ℎ: v j = Pa v i ℎ C i, i C i −1 ⋅ , ℎ . If we stack all the h ik horizontally for k = 0,…, M S − 1, we obtain the n i × n matrix ℋ(i, ⋅ ), which is i-th block row of ℋ. Intuitively, ℋ i, ⋅ is a sparse matrix with the coefficients linking the full w to w i , with zero blocks at locations whose corresponding node is v j ∈/ Pa[v i ]. The ith block-row of ℋ is of size n i × n but only has r non-empty sub-blocks, with sizes n i × n j for j ∈ {j 1 ,…, j r }, respectively. Instead, H i is a dense matrix obtained by dropping all the zero-blocks from ℋ i, ⋅ , and stores the coefficients linking w [i] to w i . The two are linked as H i w i = ℋ i, ⋅ w.
Since w = ℋw + ω, C = var w = I − ℋ −1 R I − ℋ − ⊤ , where R = b.diag{R i } and I − ℋ is block lower-triangular with unit diagonal, hence non-singular. We find the precision matrix

B.2 Properties of C −1
When not ambiguous, we use the notation X ij to denote the (i, j) block of X. An exception to this is the (i, j) block of C −1 which we denote as C −1 i, j . In SPAMTREES, C −1 i, j is nonzero if i = j or if the corresponding nodes v i and v j are connected in the moral graph G m , which is an undirected graph based on G in which an edge connects all nodes that share a child. This means that either (1) v i ∈ Pa[v j ] or vice-versa, or (2) there exists a * such that {v i , v j } ⊂ Pa[a * ]. In SPAMTREES,G m = G. In fact, suppose there is a node a * ∈ v r * such that a * ∈ Ch[v j ] ∩ Ch[v k ], where v j ∈ v r k and v k ∈ v r k . By definition of G there exists a sequence {i 1 ,…, i r * } such that Pa a* = v i 1 , ..., v i r* ⊃ v j , v k , and furthermore Pa v i ℎ = v i 1 , ..., v i ℎ − − 1 for h ≤ r * . This implies that if j = k then v j = v k , whereas if j > k then v k ∈ Pa[v j ], meaning that no additional edge is necessary to build G m .

B.2.1 EXPLICIT DERIVATION OF
, and define the "common descendants" as . Then consider a i ∈ A, v j ∈ V such that a i ∈ Pa[v j ] and denote as H i→j the matrix obtained by subsetting H j to columns corresponding to a i and note that H i j = ℋ ji . The (i, j) block of U is then if v j ∈ Ch v i Then C −1 i, j = ∑ k U ik U jk ⊤ and, as in (9), each block of the precision matrix is: where cd v i , v j = 0 implies C −1 i, j = O and I ij a zero matrix unless i = j as it is the (i, j) block of an identity matrix of dimension n × n.

B.2.2 COMPUTATION OF LARGE MATRIX INVERSES
One important aspect in building C −1 is that it requires the computation of the inverse C i −1 of dimension J i × J i for all nodes with parents, i.e. at r > 0. Unlike models which achieve scalable computations by limiting the size of the parent set (e.g. NNGPs and their blocked variant, or tessellated MGPs), this inverse is increasingly costlier when δ > 1 for Peruzzi  nodes at a higher-level of the tree as those nodes have more parents and hence larger sets of parent locations (the same conclusion holds for non-reference nodes). However, the treed structure in G allows one to avoid computing the inverse in O J i 3 . In fact, suppose we have a symmetric, positive-definite block-matrix A and we wish to compute its inverse. We write where S = C − BD −1 B ⊤ is the Schur complement of D in A. If C −1 was available, the only necessary inversion is that of S. In SPAMTREES with δ > 1, suppose v i , v j are two nodes such that Pa v j = v i ∪ Pa v i -this arises for nodes v j ∈ V r, r ≥ M δ . Regardless of whether v j is a reference node or not,

AND ITS DETERMINANT WITHOUT SPARSE CHOLESKY
Bayesian estimation of regression models requiring the computation of ZCZ ⊤ + D −1 and its determinant use the Sherman-Morrison-Woodbury matrix identity to find Cholesky factorization of C −1 + ∑ can be used as typically Σ is diagonal or block-diagonal, thus maintaining the sparsity structure of C −1 . Sparse Cholesky libraries (e.g. Cholmod, Chen et al., 2008), which are embedded in software or high-level languages such as MATLAB™ or the Matrix package for R, scale to large sparse matrices but are either too flexible or too restrictive in our use cases: (1) we know G and its properties in advance; (2) SPAMTREES take advantage of block structures and grouped data. In fact, sparse matrix libraries typically are agnostic of G and heuristically attempt to infer a sparse G given its moralized counterpart. While this operation is typically performed once, a priori knowledge of G implies that reliance on such libraries is in principle unnecessary.
We thus take advantage of the known structure in G to derive direct algorithms for computing C −1 + ∑ −1 and its determinant. In the discussion below we consider δ = M , Peruzzi  noting here that choosing δ = 1 simplifies the treatment as cd(v i , v j ) = v i if v i = v j , and it is empty otherwise. We now show how (21) leads to Algorithm 4 for the decomposition of any precision matrix Λ which conforms to G -i.e. it has the same block-sparse structure as a precision matrix built as in Section B.1. Suppose from Λ we seek a block lower-triangular matrix L and a block diagonal D such that Start with v i , v j taken from the leaf nodes, i.e. v i , v j ∈ V M . Then cd v i , v j = 0 and we set we then set L ii = O and get the i-th block of D i simply setting D i = Λ ii . Proceeding where L ij and D i have been fixed at the previous step; this results in D j = Λ jj − L ij Then, the s-th (of M ) step takes v j ∈ V M − s ∩ Pa v i and v i ∈ V M − s + 1, , implying cd v i , v j = v i ∪ ch v i . Nothing that F * = ∑ v k ∈ ch v i I ki − L ki ⊤ D k I kj − L kj has been fixed at previous steps since each v k is at level M−s + 2, we split the sum in (21)  In practice, a block matrix with K 2 blocks can be represented as a K 2 array with rows and columns indexed by nodes in G and matrix elements which may be zero-dimensional whenever corresponding to blocks of zeros. The specification of all algorithms in block notation allows us to never deal with large (sparse) matrices in practice but only with small block matrices indexed by nodes in G, bypassing the need for external sparse matrix libraries. Specifically we use the above algorithms to compute Λ −1 = C −1 + Σ −1 and its determinant: Λ −1 = I − L −1 D −1 I − L − ⊤ and Λ −1 = i = 1 M S 1/ D i . We have not distinguished non-reference and reference nodes in this discussion. In cases in which the non-reference set is large, we note that the conditional independence of all non-reference locations, given their parents, results in C −1 i, i being diagonal for all ℓ ∈ U (i.e. η(ℓ) = v i ∈ B). This portion of the precision matrix can just be stored as a column vector.

B.2.4 SPARSITY OF C −1
We calculate the sparsity in the precision matrix; considering an enumeration of nodes by level in G, denote n ij = η −1 v ij , m j = V j , and J ij = η −1 Pa v ij , and noting that by We outline recursive properties of C induced by G when δ = M . In the case 1 < δ < M , these properties hold for nodes at or above level M δ , using A M δ as root. We focus on paths in G. These can be represented as sequences of nodes v i 1 , …, v i r such that v i j , …, v i k ⊂ P a v i k + 1 for 1 < j < k < r. Take two successive elements of such a sequence, Consider E w j | w [j] = H j w j = C j, j C j −1 w j ET and R j = var w j | w j = C j, j − C j, j C j −1 C j , j . By (22) we can write Now define the covariance function K i ℓ, ℓ′ = C ℓ, ℓ′ − C ℓ, i C i −1 C i , ℓ′ ; recalling that the reference set is S = ∪ i = 0 M − 1 ∪ j = 1 m j S j we use a shorthand notation for these subsets: the conditionals governing the process as prescribed by G. We can also write the above as E w j | w j = s = i 0 i r cov e ℎ , e k = E cov e ℎ , e k w ℎ , w ℎ + cov E e ℎ w ℎ , w ℎ = cov E e ℎ w ℎ , w ℎ , E e k w ℎ , w ℎ = 0.
The above results also imply C j, j C j −1 C j , j = s = i 0 i r K s j, s K s −1 s, s Ks s, j and suggest an additive representation via orthogonal basis functions: Finally, considering the same sequence of nodes, recursively introduce the covariance functions F 0 (r, s) = C r,s and for j > 1, F j r, s = F j − 1 r, s − F j − 1 r, j − 1 F j − 1 −1 j − 1, j − 1 F j − 1 j − 1, s .
We get F j + 1 r, s = F j r, s − F j r, j F j −1 j, j F j j, s using 22 = F j − 1 r, s − F j − 1 r, j − 1: j F j − 1 −1 j − 1: j , j − 1: j F j − 1 −1 j − 1: j , s = C r, s − C r, 0: j C −1 0: j , 0: j , s = C r, s − C r, j + 1 C j + 1 −1 C j + 1 , s which can be iterated forward and results in an additional recursive way to compute covariances in SPAMTREES. Notice that while K j is formulated using the inverse of J j × J j matrix C [j] , the F j 's require inversion of smaller n j × n j matrices F j−1 (j-1, j-1).

B.4.1 δ = 1
Choosing depth δ = 1 results in each node having exactly 1 parent. In this case the path P k j = v i 1 , …, v i r from v k to v j , where v i 1 = v k , v ir = v j and v i ℎ = P a v i ℎ + 1 , is unique, and there is thus no distinction between shortest and longest paths: P k j = P k j = P k j . Then denote H ⃛ k j = H i r ⋅ H i r − 1 ⋯H i 1 . Let v z be the concestor between v i and v j i.e. v z = con v i , v j = argmax v k ∈ V k: P k i ∩ P k i ≠ 0 and the associated paths P z i = v i 1 , …, v i r i and P z j = v j 1 , …, v j r j where v i 1 = v j 1 = v z , v i r i = v i and v j r j = v j . Then we can write where v i r i ∼ N 0, R i r i and proceed expanding w i r i − 1 to get w i r i = H i r i H i r i − 1 w i r i − 2 + v i r i − 1 + v i r i = H i r i H i r i − 1 w i r i − 2 + H i r i v i r i − 1 + v i r i w i = H i r i ⋯H i 1 w i 1 + v i = H ⃛ z i w z + v i where v i is independent of w z . After proceeding analogously with w j , take ℓ i , ℓ j such that η(ℓ i ) = v i and η(ℓ j ) = v j . Then Cov p w ℓ i , w ℓ j = H ⃛ z i ℓ i C z H ⃛ z j ℓ j ⊤ , (24) where H ⃛ z i ℓ i = C ℓ i , S i C i −1 H ⃛ z i and similarly for H ⃛ z j ℓ j .

B.4.2. 1 < δ < M
Take two nodes v i ,v j ∈ V . If Pa v i ∩ Pa v j ≠ 0 then we apply the same logic as in B.4.3 using v z = con(v i , v j ) as root.
where v x ∈ A M δ is the parent node of v i at level M δ . The final result of (14) is then achieved by noting that the relevant subgraph linking v x and v j has depth δ x = 1 and thus Cov(w x , w j ) can be found via B.4.1, then Cov w i , w j = C ix C x −1 Cov w x , w y = F i Cov w x , w y . Notice that F i directly uses the directed edge v x v i in G; for this reason the path between w i and w i = con w x , w j is the actually the shortest path and we have v z ⋯ v x v i .

B.4.3 δ = M
Take v i , v j ∈ V and the full paths from the root P 0 i = i 0 , …, i r i and P 0 j = j 0 , …, j r j , respectively. Then using (23)  where Cov e i , e j = 0. Then since e s are independent and e s N 0, K s s, s we find Cov p w i , w j = Cov s ∈ P 0 i ∩ P 0 j K s i, s K s −1 s, s e s + e i , s ∈ P 0 i ∩ P 0 j K s j, s K s −1 s, s e s + e j = s ∈ P 0 i ∩ P 0 j K s i, s K s −1 s, s K s s, j + 1 i = j K i i, i .
We conclude by noting that δ = M implies P 0 i ∩ P 0 j = Pa v i ∩ Pa v j ; considering two locations ℓ i , ℓ j ∈ D * such that η ℓ i = v i and η ℓ j = v j we obtain Cov p w ℓ i , w ℓ j = s: v s ∈ Pa v i ∩ Pa v j K s ℓ i , s K s −1 s, s K s s, ℓ j + 1 ℓ i = ℓ j K i ℓ i , ℓ j .

B.5. Computational cost
We make some assumptions here to simplify the calculation of overall cost: first, we assume that reference locations are all observed S ⊂ T, and consequently U = T\S. Second, we assume that all reference subsets have the same size i.e.   The storage requirements are driven by the covariance at parent locations C j for nodes v j with Pa v j ≠ 0 i.e. all reference nodes at level r = 1, …, M − 1 and non-reference nodes.
Taking δ = M , suppose v i is the last parent of v j , meaning v i ∪ Pa v i = Pa v j . Then proceeding analogously as above. This stepwise procedure is stopped when either the tree reaches a predetermined height M , or when there is an insufficient number of remaining locations to build reference subsets of size n S . The remaining locations are assigned to the leaf nodes via η B as defined in Section 2.1 in order to include at least one neighboring realization of the process from the same variable.
One specific issue arises when multivariate data are imbalanced, i.e. one of the margins is observed at a much sparser grid, e.g. in Section 4.2 PRCP is collected at a ratio of 1:10 locations compared to other variables. In these cases, if locations were chosen uniformly at random to build the reference subsets then the root nodes would be associated via η to reference subsets which likely do not contain such sparsely observed variables. This scenario goes against the intuition of 1 suggesting that a naïve approach would result in poor performance at the sparsely-observed margins. To avoid such a scenario, we bias the sampling of locations to favor those at which the sparsely-observed variables are recorded. As a result, in Section 4.2 near-root nodes are associated to reference subsets in which all variables are balanced; the imbalances of the data are reflected by imbalanced leaf nodes instead.
The source code for SPAMTREES is available at https://CRAN.R-project.org/ package=SpamTree and can be installed as an R package. The SPAMTREE package is written in C++ using the Armadillo library for linear algebra (Sanderson and Curtin, 2016) interfaced to R via RcppArmadillo (Eddelbuettel and Sanderson, 2014). All matrix operations are performed efficiently by linkage to the LAPACK and BLAS libraries (Blackford et al., 2002;Anderson et al., 1999) as implemented in OpenBLAS 0.3.10 (Zhang, 2020) or the Intel Math Kernel Library. Multithreaded operations proceed via OpenMP (Dagum and Menon, 1998).

C.1. On the dependencies on sparse Cholesky libraries
SPAMTREES do not require the use of external Cholesky libraries because Cholesky-like algorithms can be written explicitly by referring to the treed DAG and its edges. This is unusual for DAG-based models commonly used in geostatistical settings. For example, an NNGP model uses neighbors to build a DAG. The DAG can be used to fill the L and D matrices leading to ℓDℓ ⊤ = C −1 , where C −1 is the sparse precision matrix of the latent process. When one then adds measurement error in a regression setting and marginalizes out the latent process, the goal is to find the Cholesky decomposition of C −1 + τ 2 I n . The original NNGP DAG used for L and D is not useful for this purpose -hence the need for NNGP to use sparse Cholesky libraries in collapsed samplers . On the other hand, with SPAMTREES we can still look at the original DAG to "update" Land D by using Algorithm 4. We included it in the Appendix as our software package implements the Gibbs sampler in the main article, which does not involve Cholesky decompositions of large sparse precision matrices.
Furthermore, one of the initial steps in Cholesky algorithms for sparse symmetric positivedefinite matrices involves finding "good" reordering rows and columns. These reorderings simplify the (undirected) graphical model that corresponds to the sparsity structure in the matrix. Once a simple-enough graphical model is found heuristically, it is used for the decomposition. On the other hand, the sparse Cholesky algorithm for SPAMTREES can be written explicitly using the underlying DAG, without any intermediate step, because it is fixed and with a convenient treed structure. It might be possible to write software that out-performs excellent libraries such as CHOLMOD (Chen et al., 2008) at decomposing the matrices needed in collapsed sampling algorithms for SPAMTREES.
Regarding a more general perspective about dependencies on well established software libraries, we should clarify that the software for SPAMTREE does take advantage of highly efficient libraries such as BLAS/LAPACK provided in Intel MKL 2019.5, OpenMP (Dagum and Menon, 1998) for parallelizing the algorithms as described in the main article. These libraries are optional and our code does not strictly depend on them. For example, noting that it is considerably more difficult to compile OpenMP code on Macs, one can just disable OpenMP and let the Accelerate BLAS (rather than OpenBLAS or Intel MKL) deal with all matrix algebra for Apple computers, at the cost of some performance in big data settings.
Our code also has R package dependencies for data pre-processing, but these are peripheral to the proposed methods and algorithms. Any improvement in the upstream libraries we used for coding SPAMTREE will positively impact the performance of our software.

C.2.1. SIMULATED DATASETS
SPAMTREES with full depth are implemented by targeting reference subsets of size n S = 25 and tress with c = 4 additional children for each branch. The tree is built starting from a 2 × 2 partition of the domain, hence there are 4 root nodes with no parents in the DAG. The cherry-pickying function η is set as in Section 2.1; with these settings the tree height is M = 3. For SPAMTREES with depth δ = 1 we build the tree with reference subsets of size n S = 80 and c = 4. Multivariate Q-MGPs are implemented via axis-parallel partitioning using 57 intervals along each axis. The multivariate SPDE-INLA method was implemented following the examples in Krainski et al. (2019), Chapter 3, setting the grid size to 15 × 15 to limit the compute time to 15 seconds when using 10 CPU threads. BART was implemented on each dataset via the wbart function in the R package BART; the set of covariates for BART was built using the spatial coordinates in addition to a binary variable representing the output variable index (i.e. taking value 1 whenever y i is of the first outcome variable, 0 otherwise).

C.2.2. MODIS-TERRA AND GHCN
The implemented SPAMTREES are built with 36 root nodes and c = 6 additional children for each level of the tree, 25 reference locations for each tree node, and for up to M = 5 levels of the tree and δ = 5 (i.e. full depth). The non-reference observed locations are linked to leaves via cherry-pickying as in Section 2.1. Multivariate models were run on an AMD Epyc 7452based virtual machine with 256GB of memory in the Microsoft Azure cloud; the SPAMTREE R package was set to run on 20 CPU threads, on R version 4.0.3 linked to the Intel Math Kernel Library (MKL) version 2019.5-075. The univariate models were run on an AMD Ryzen 5950X-based dedicated server with 128GB of memory, on 16 threads, R version 4.1.1 linked to Intel MKL 2019.5-075. The univariate NNGP model was implemented using R package spNNGP  using 20 neighbors for all outcomes and a "latent" algorithm. The univariate SPAMTREE was implemented on each outcome with full depth, c = 16 additional chidren for each level of the tree, and 25 reference locations for each node.
The MGP model was implemented via the development package at github.com/mkln/ meshgp targeting a block size with 4 spatial locations, resulting in an effective average block dimension of 20. Caching was unavailable due to the irregularly spaced PRCP values. Fewer MCMC iterations were run compared to SPAMTREES to limit total runtime to less than 16h. Prediction and estimation performance on multivariate synthetic data. The four columns on the right refer to root mean square error (RMDE) and mean absolute error (MAE) in out-of-sample predictions, average coverage of empirical 95% prediction intervals, and RMDE in the estimation of θ.